Houdiniメモ : 1次元の波動方程式を実装する

Solverノードを使って1次元の波動方程式を実装してみました。

https://gyazo.com/847c475ec84a4e262df711834fb8db4c

ネットワーク解説

■Lineノード

LineノードではX方向のラインを作成します。

https://gyazo.com/90813142a26139b50486f8cc323d34ac

■Resampleノード

https://gyazo.com/4c6272cd803d879fafe70f5c1573f350

■Transformノード

ラインの適当なポイントを持ち上げてます。 この状態が波の初期状態となります。

https://gyazo.com/7b2785ec353ff3afd756bf94d8454ec9

Solverノードの内部実装

Solverノードの内部ではAttributeWrangleノードを二つ作成します。

■1つ目のAttributeWrangleノード(速度の計算)

Prev_FrameノードにAttributeWrangleノードを接続し、速度を計算するVEXを記述します。

https://gyazo.com/368c50918ef53761a68054cdfccbff55

code:calc_speed(c)

// 波の端は動かさない(端点は固定)

if (@ptnum < 1 || @ptnum >= @numpt - 1) { return; }

vector p0 = point(geoself(), "P", @ptnum - 1);

vector p1 = point(geoself(), "P", @ptnum);

vector p2 = point(geoself(), "P", @ptnum + 1);

// 時刻tにおける座標xにあるポイントのY座標を波の変異 u(x, t)とする

float u_prev = p0.y; // u(x - dx, t)

float u_center = p1.y; // u(x, t)

float u_next = p2.y; // u(x + dx, t)

// 波の加速度の計算

float dx = (p2.x - p1.x); // xの微分

float dt = 1.0 / (float)$FPS; // 時間の微分

float v = chf("v"); // 波が伝わる速さ

f@Accel = v * v * (u_next + u_prev - 2 * u_center) / (dx * dx); // 微分の近似計算で加速度を求める

// 加速度を速度に加える

f@Speed += f@Accel * dt;

■2つ目のAttributeWrangleノード(ポイントの座標更新)

下流のAttributeWrangleでは、計算した速度を元にポイントを動かします。

code:move(c)

float dt = 1.0 / (float)$FPS;

@P.y += f@Speed * dt;

https://gyazo.com/6556bdf6acd50d9fba354e275d8ae520

参考

波動方程式の数値シミュレーション(筑波大学)